(*Export["GeoNu.pdf",%]*)

(*
Solar neutrino spectrum data from http://www.sns.ias.edu/~jnb/SNdata/sndata.html 
Fluxes from Table 2, last row, http://arxiv.org/abs/astro-ph/0412440
Two be7 lines intensities are 10% and 90%.
*)

ppSRC=(Drop[Import["http://www.sns.ias.edu/~jnb/SNdata/Export/PPenergyspectrum/ppenergytab","Table"],3]//Transpose)/.{0.`->10.0^-12};

ppDAT=Transpose[{Table[ppSRC[[ii]],{ii,1,7,2}]//Flatten,Table[ppSRC[[ii]],{ii,2,8,2}]//Flatten}];

b8DAT=Drop[Drop[Import["http://www.sns.ias.edu/~jnb/SNdata/Export/B8spectrum/b8spectrum.txt","Table"],16],-3];

b8=Table[{b8DAT[[ii,1]],4.59*10^6*b8DAT[[ii,2]]},{ii,1,Dimensions[b8DAT][[1]]}] ;

pp=Table[{ppDAT[[ii,1]],6.05*10^10*ppDAT[[ii,2]]},{ii,1,Dimensions[ppDAT][[1]]}];

n13DAT=Import["http://www.sns.ias.edu/~jnb/SNdata/Export/CNOspectra/n13.dat","Table"]/.{0.`->10.0^-12};
o15DAT=Import["http://www.sns.ias.edu/~jnb/SNdata/Export/CNOspectra/o15.dat","Table"]/.{0.`->10.0^-12};
f17DAT=Import["http://www.sns.ias.edu/~jnb/SNdata/Export/CNOspectra/f17.dat","Table"]/.{0.`->10.0^-12};

hepDAT=Import["http://www.sns.ias.edu/~jnb/SNdata/Export/Hepspectrum/hepspectrum.dat","Table"];

n13={Transpose[n13DAT][[1]],Transpose[n13DAT][[2]]*2.3*10^8}//Transpose;
o15={Transpose[o15DAT][[1]],Transpose[o15DAT][[2]]*1.47*10^8}//Transpose;
f17={Transpose[f17DAT][[1]],Transpose[f17DAT][[2]]*3.31*10^6}//Transpose;

hep={Transpose[hepDAT][[1]],Transpose[hepDAT][[2]]*8.23*10^3}//Transpose;

be7SRC=Import["http://www.sns.ias.edu/~jnb/SNdata/Export/7Belineshape/7belineshape.dat"]/.{0.`->10.0^-12};

be7DATa=Drop[Drop[be7SRC,3],-87];

Dimensions[be7DATa]

be7DATb=Drop[be7SRC,99];

Dimensions[be7DATb]

be7a=Append[Table[{(be7DATa[[ii,1]]+861.3)/1000,0.9*4.38*10^9*be7DATa[[ii,2]]},{ii,1,92}],{0.86873,0.1}];

be7b=Append[Table[{(be7DATb[[ii,1]]+384.3)/1000,0.1*4.38*10^9*be7DATb[[ii,2]]},{ii,1,83}],{0.39169+$MachineEpsilon,0.1}];

solar=ListLogLogPlot[{pp,b8,n13,o15,f17,be7a,be7b,hep},Joined->True,PlotRange->{{0.099,20},{10^1,10^12}},Frame->True,PlotStyle->{Red,Blue,{Dashed,Red},{Dashed,Green},{Dashed,Blue},Black,Black,Yellow},BaseStyle->FontSize->16,FrameLabel->{"(Subscript[E, ν])[MeV]","(Subscript[F, ν])[(cm^(-2))(s^(-1))(MeV^(-1))]"},PlotLegends->{"pp","\!\(\*SuperscriptBox[\(\\ \), \(8\)]\)B","CNO\!\(\*SuperscriptBox[\(\\ \), \(13\)]\)N","CNO\!\(\*SuperscriptBox[\(\\ \), \(15\)]\)O","CNO\!\(\*SuperscriptBox[\(\\ \), \(17\)]\)F","\!\(\*SuperscriptBox[\(\\ \), \(7\)]\)Be","\!\(\*SuperscriptBox[\(\\ \), \(7\)]\)Be","hep"},ImageSize->600]

(*Export["SolarNu.pdf",%]*)

Q=2ParticleData["Proton","Mass"]-IsotopeData["Hydrogen2","AtomicMass"]//QuantityMagnitude

IsotopeData["Hydrogen2","BindingEnergy"]

UnitConvert[2IsotopeData["Hydrogen1","AtomicMass"]-IsotopeData["Hydrogen2","AtomicMass"],"MeV/c^2"]

Select[ParticleData["Baryon"],ParticleData[#,"Stable"]==True&]

2ParticleData["Proton","Mass"]-ParticleData[,"Mass"]

me=ParticleData["Electron","Mass"]//QuantityMagnitude

Show[ListPlot[pp],Plot[1.17*10^13 Eν^2(Q-Eν+me)Sqrt[(Q-Eν+me)^2-me^2],{Eν,0,Q}]]

(*Export["SolarNu.pdf",%]*)

(* 
Geoneutrino Spectrum and Luminosity from Sanhiro Enomoto PhD.
http://www.awa.tohoku.ac.jp/~sanshiro/research/geoneutrino/spectrum/index.html
*)

geonuU238=Drop[Import["http://www.awa.tohoku.ac.jp/~sanshiro/research/geoneutrino/spectrum/AntineutrinoSpectrum_238U.knt","Table"],6];
geonuTh232=Drop[Import["http://www.awa.tohoku.ac.jp/~sanshiro/research/geoneutrino/spectrum/AntineutrinoSpectrum_232Th.knt","Table"],6];
geonuK40=Drop[Import["http://www.awa.tohoku.ac.jp/~sanshiro/research/geoneutrino/spectrum/AntineutrinoSpectrum_40K.knt","Table"],6];

NormK40=Sum[geonuK40[[ii]][[2]],{ii,1,Dimensions[geonuK40][[1]]}]
NormU238=Sum[geonuU238[[ii]][[2]],{ii,1,Dimensions[geonuU238][[1]]}]
NormTh232=Sum[geonuTh232[[ii]][[2]],{ii,1,Dimensions[geonuTh232][[1]]}]

fluxU238=2.34*10^6 (* fluxes at Kamioka/Japan *)
fluxTh232=1.99*10^6

fluxK40=0.5*10^6 (* dla ilustracji, wartość wzięta z powietrza!*)

U238=
Table[ {geonuU238[[ii]][[1]] /1000,fluxU238*geonuU238[[ii]][[2]]/NormU238*1000},{ii,1,Dimensions[geonuU238][[1]]}];

Th232=
Table[ {geonuTh232[[ii]][[1]] /1000,fluxTh232*geonuTh232[[ii]][[2]]/NormTh232*1000},{ii,1,Dimensions[geonuTh232][[1]]}];

K40=
Table[ {geonuK40[[ii]][[1]] /1000,fluxK40*geonuK40[[ii]][[2]]/NormK40*1000},{ii,1,Dimensions[geonuK40][[1]]}];

geonu=Table[ {geonuTh232[[ii]][[1]] /1000,fluxTh232*geonuTh232[[ii]][[2]]/NormTh232*1000+fluxU238*geonuU238[[ii]][[2]]/NormU238*1000},{ii,1,Dimensions[geonuTh232][[1]]}];

dsnb=LogLogPlot[10Eν^2/(1+Exp[Eν/(kT)])/(1.8030853547393915` (kT^3))/.kT->4,{Eν,0.1,100},PlotRange->{0.1,10^12},PlotStyle->{Black,DotDashed},PlotLegends->{"DSNB"}]

geo=Show[ListLogLogPlot[{U238,Th232,K40},PlotRange->{{0.0999,20.0},{10^-1,10^9}},Joined->True,PlotStyle->{Red,Blue,Green},PlotLegends->{"\!\(\*SuperscriptBox[\(\), \(238\)]\)U","\!\(\*SuperscriptBox[\(\), \(232\)]\)Th","\!\(\*SuperscriptBox[\(\), \(40\)]\)K"}],dsnb,AxesLabel->{"(Subscript[E, ν])[MeV]","(Subscript[F, ν])[(cm^(-2))(s^(-1))(MeV^(-1))]"},BaseStyle->FontSize->16,ImageSize->600]